Chiral magnetism and helimagnons in a pyrochlore antiferromagnet 
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Recent neutron scattering measurements on the spinel CdCr2 04 revealed a rare example of helical 
magnetic order in geometrically frustrated pyrochlore antiferromagnet. The spin spiral character- 
ized by an incommensurate wavevector Q = 2n(0, S, 1) with 8 ~ 0.09 is accompanied by a tetragonal 
distortion. Here we conduct a systematic study on the magnetic ground state resulting from the 
interplay between the Dzyaloshinskii-Moriya interaction and further neighbor exchange couplings, 
two of the most important mechanisms for stabilizing incommensurate spin orders. We compute 
the low-energy spin-wave spectrum based on a microscopic spin Hamiltonian and find a dispersion 
relation characteristic of the helimagnons. By numerically integrating the Landau-Lifshitz-Gilbert 
equation with realistic model parameters, an overall agreement between experiment and the nu- 
merical spectrum, lending further support to the view that a softened optical phonon triggers the 
magnetic transition and endows the lattice a chirality. 



I. INTRODUCTION 

Magnets with geometrical frustration provide a fertile 
ground for studying complex spin orders and unconven- 
tional magnetic phases. 1 The defining feature of strong 
geometrical frustration is the occurrence of extensively 
degenerate ground states. This usually arises when the 
arrangement of spins on a lattice precludes satisfying all 
interactions at the same time. A great deal of attention 
has been devoted to the so-called pyrochlore magnets in 
which the magnetic ions form a three-dimensional con- 
nected network of corner-sharing tetrahedra. 2 With mag- 
netic interactions restricted to only nearest neighbors, 
classical spins on the pyrochlore lattice remain disordered 
even at temperatures well below the Curie- Weiss energy 
scale. Magnetic fluctuations in this liquid-like phase are 
subject to strong local constraints that maintain van- 
ishing total spin on each tetrahedron. This disordered 
yet highly correlated phase is often called a cooperative 
paramagnet. Because of the huge degeneracy, the magnet 
is sensitive to nominally small perturbations. The low- 
temperature ordering thus depends subtly on the residual 
perturbations present in the system. 

The recent discovery of a novel chiral spin structure in 
spinel CdCr2C>4 has generated much interest both theo- 
retically and experimentally. 3,4 The magnetic Cr ions in 
this compound have spin S — 3/2 and no orbital degrees 
of freedom, providing a very good realization of Hciscn- 
bcrg antiferromagnet on the pyrochlore lattice. Despite 
a rather high Curie- Weiss temperature |©cw| ~ 70 K, 
a spiral magnetic order sets in only at Xjv = 7.8 K, in- 
dicating a high degree of frustration. The discontinuous 
magnetic transition at TV is accompanied by a struc- 
tural distortion which lowers the crystal symmetry from 
cubic to tetragonal (a = b < c) with an elongated unit 
cell: (c — a)/c w 5 x 10~ 3 . Polarized neutron-scattering 
experiment reveals an incommensurate magnetic order 
in which coplanar spins rotate about either a or b axes 
with a period of roughly ten lattice constants; the or- 
dering wavevector Q = 2tt(0,S, 1) with S « 0.09. It is 



worth noting that a very similar coplanar helical order 
has also been reported in the weak itinerant antiferro- 
magnet YMn 2 5,6 . 

Spin-lattice coupling has been shown to play an 
important role in relieving magnetic frustrations in 
several chromium spinels. 7-10 In particular, there are 
strong experimental evidences indicating that the mag- 
netic transition in CdCr204 is mainly driven by spin- 
lattice coupling. These include anomalies of the clastic 
constants, 11,12 splitting of optical phonons, 13,14 and en- 
ergy shift of zone-boundary phonons 15 at the magnetic 
ordering temperature. The spin-lattice or magnctoclastic 
coupling arises from the dependence of exchange on the 
ion displacements fry- of spins: (dJ/dr)(Si ■ Sj) fry. 16,17 
As a result, the lattice distortion creates disparities be- 
tween nearest-neighbor exchange constants and paves the 
way for magnetic ordering. However, spin-lattice cou- 
pling usually favors collincar magnetic orders. Indeed, 
assuming a simple elastic energy cost k(Srij) 2 /2 and in- 
tegrating out the displacements 5rtj in the above model 
gives rise to a biquadratic spin interaction: —(Si • Sj) 2 , 
which is minimized by collinear spins. While the mag- 
netic frustration is relieved by the distortion, the appear- 
ance of helical spin order in CdCr204 thus comes from 
other mechanisms. 

Magnetic spirals resulting from the competition be- 
tween nearest neighbor and next nearest neighbor in- 
teractions are a common theme in frustrated spin 
systems. 18,19 However, this is not the case for pyrochlore 
lattice. Detailed mean-field and numerical investiga- 
tions show that inclusion of further-neighbor interactions 
partially relieves the frustration and sometimes selects 
rather complicated multiple-q spin orders, but fails to 
stabilize magnetic spirals with a definite handedness. 20,21 
Another mechanism for the formation of spirals is the 
Dzyaloshinskii-Moriya (DM) interaction which originates 
from the relativistic spin-orbit coupling. 22,23 Although 
symmetry considerations allow the presence of DM terms 
in the pyrochlore lattice, 24 the ground states selected by 
the DM interaction are non-chiral magnetic orders with 



wavcvector q = 0. 

The origin of the spiral order as well as the lattice dis- 
tortion in CdCi'204 has been consistently explained in 
Ref. 27. In this scenario, the magnetic transition is trig- 
gered by the softening of q = optical phonons with odd 
parity. The two types of tetrahedra with different orien- 
tations are flattened along the a and b directions, respec- 
tively, resulting in an overall tetragonal elongation along 
the c axis. The lattice distortion relieves the magnetic 
frustration by creating disparities among the nearest- 
neighbor exchange constants and stabilizes a collinear 
Neel order with wavevector Q = 2n(0, 0, 1). Moreover, 
the broken parity due to the staggered distortion gives 
rise to a chiral pyrochlore structure. The collinear spins 
are twisted into a long-period spiral with Q = 27r(0, 6, 1) 
as the structural chirality is transferred to magnetic or- 
der by relativistic spin-orbit interaction. 17,27 The broken 
parity in CdCr204 has also been confirmed by recent op- 
tical measurements. 14 

In this paper we undertake a systematic study of the 
general helical order resulting from the interplay between 
DM interaction and further-neighbor couplings in py- 
rochlore antifcrromagncts. We also employ various nu- 
merical methods to investigate the spin-wave excitations 
in the spiral magnetic orders. We first revisit the contin- 
uum theory which is valid in the J — > 00 limit and show 
that the resultant Hamiltonian possesses a continuous 
symmetry similar to that in smectic liquid crystals. 2 
In particular, despite the broken symmetry of the heli- 
cal order is described by a 0(3) order parameter, the 
three associated Goldstone modes reduce to a single 
smectic-like phonon mode which emerges on scales longer 
than the helical pitch. This low-energy gapless excita- 
tions, dubbed helimagnons in Ref. 31, exhibit a highly 
anisotropic dispersion relation. 

The emergent continuous symmetry of the effective 
theory is broken by the anisotropic crystal fields at the 
microscopic level. We first adopt the spherical approx- 
imation to investigate the resultant magnetic ground 
states. Three typical helical orders with wavevectors 
(0, 0, 1 — (5), (0, S, 1) and (6, S, 1) are stabilized depending 
on the relative strengths of DM interaction D and third 
neighbor interaction J3. We then focus on the (0,6,1) 
spiral, which is relevant to the case of CdCi"2 04, and con- 
sider perturbation corrections to the coplanar spiral solu- 
tion obtained from the continuum theory. We show that 
the first-order corrections do not shift the spiral pitch but 
introduces a small out-of-plane spin component. 

Next we study the magnetic excitations in the spi- 
ral magnetic ground state. The incommensurate nature 
of the helical order presents a challenging task for the 
calculation of spin-wave spectrum. To circumvent this 
problem we employ the conventional Holstein-Primakoff 
method with a large commensurate unit cell. We also 
perform dynamics simulations based on the Landau- 
Lifshitz-Gilbert (LLG) equations. 32 The LLG simula- 
tions with a large damping also provide an efficient ap- 
proach to obtain the accurate spin structures in the equi- 
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FIG. 1: (Color online) Schematic picture of (a) even and (b) 
odd tetragonal distortions with E g and E u symmetries, re- 
spectively. Both types of distortion preserve the lattice trans- 
lational symmetry and could result in an overall elongation 
along the c axis. Panel (c) shows the DM vectors Dij on the 
six inequivalent NN bonds of the pyrochlore lattice. 

librium state. The magnon spectrum obtained from our 
dynamics simulation exhibits a gapless mode at the heli- 
cal wavevector Q = 27r(0, 6, 1) reminiscent of the smectic- 
like phonon mode or helimagnons mentioned above. ' 
An overall agreement is obtained between the numerical 
calculations and the experimental data on CdCr2 04. 

The remainder of this paper is organized as follows. 
In Sec. II we consider the helical magnetic order in py- 
rochlore antifcrromagnct. After presenting a minimal 
model relevant for CdCr204, we first review the contin- 
uum theory of the magnetic spirals on pyrochlore. We 
then demonstrate the existence of a continuous symme- 
try similar to that in smectic liquid crystals. We obtain 
a phase diagram of the helical order based on the spher- 
ical approximation. Sec. Ill is devoted to a discussion of 
the magnon excitations of pyrochlore magnet. A system- 
atic analysis of the low-energy acoustic mode indicates 
a highly anisotropic dispersion characteristic of that of 
helimagnons. We then perform dynamics simulation to 
obtain the spin-wave spectrum of the helical spin order 
with wavevector Q = 2n(0, 5, 1) and compare the results 
with available experimental data. We conclude in Sec. V 
with a summary and discussion of our results. 



II. HELICAL MAGNETIC ORDER 

Our starting point is a classical Hciscnbcrg model on 
the pyrochlore lattice described by the Hamiltonian 

H = ( J + K n) s * • S, • X D ' iS - x s j) 

■ ./•> X s < • s < • E s <- s <- 

<«» <««») 

The dominant term is the nearest-neighbor (NN) ex- 
change interaction with an antiferromagnctic sign J > 0. 
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Due to the special geometry of the pyrochlore lattice, 
this term can be recast to ( J/2) |M(r)| 2 , where M 
denotes the vector sum of spins on a tetrahedron at r. 
Minimization of this term requires the vanishing of to- 
tal magnetic moment on every tetrahedron but leaves an 
extensively degenerate manifold. 

The exchange anisotropy Kij results from a distorted 
lattice through magnetoelastic coupling. The structural 
distortion in CdCr204 reduces the crystal symmetry 
from cubic to tetragonal while preserving the lattice 
translational symmetry. 3 Here we consider q = lat- 
tice distortions with even and odd parities; see Fig. 1(a) 
and (b). The dominant E u component results in a stag- 
gered distortion: tetrahedra of types I and II are flattened 
along the a and b directions, respectively; the lattice is 
elongated overall, a — b < c. The resulting space group, 
14x22, lacks the inversion symmetry. Fig. 2 shows the 
exchange anisotropy of the six NN bonds in the unit cell. 
As pointed out in Ref. 16, this tetragonal distortion com- 
pletely relieves the magnetic frustration and stabilizes 
a collinear Neel order with wavevector Qo = 27r(0, 0,1) 
shown in Fig. 2. 

The DM term • (Si x Sj) originates from the rel- 
ativistic spin-orbit interaction A(L ■ S). This term is al- 
lowed on the ideal pyrochlore lattice, where the bonds are 
not centrosymmctric as required by the so-called Moriya 
rules. 23 Moreover, the high symmetry of the pyrochlore 
lattice completely determines the orientations of the DM 
vectors [see Fig. 1(c)] up to a multiplicative factor. Fi- 
nally, we also include second and third exchanges cou- 
plings J2 and J3, respectively, in the model Hamiltonian. 
First principle ab initio calculations indicate that these 
further-neighbor exchanges are important for the model- 
ing of spinel CdCr2 04. 27,34 In the limit dominated by a 
large NN exchange, the ground-state properties actually 
depend only on the difference between J2 and J3. The 
relative shift in energy for any pair of ground states due 
to a small J3 is identical to the effect of a J2 of the same 
magnitude with an opposite sign. 21 



A. Gradient approximation and emergent 
continuous symmetry 

We first briefly review the continuum theory of helical 
orders derived from the Qo = 27r(0, 0, 1) Neel order. The 
effective energy functional is obtained based on a gradi- 
ent approximation of the staggered magnetizations. 17,27 
The antifcrromagnetic order on the pyrochlore lattice is 
characterized by three order parameters: 



Lj = (So + S x - S 2 - S 3 )/4, 
L 2 = (So - Sx + S 2 - S 3 )/4, 
L 3 = (So - Sx - S 2 + S 3 )/4, 

where S p denotes magnetic moments at /xth sublattice 
For example, the collinear Neel state shown in Fig. 2 is 
described by order parameters: M = 0, L2 = L3 = 0, 
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FIG. 2: (Color online) A Neel order on the pyrochlore lat- 
tice with wavevector q = 2n(0, 0, 1). This collinear spin or- 
der is stabilized by a staggered distortion with E u symmetry. 
Tetrahedra of type I and II are flattened along the a and b 
directions, respectively. The exchange constants on various 
NN bonds are J + K u - K g (solid), J - 2K U ~ K g (dashed), 
and J + K u + 2K g (dotted). 



and Lx = ni e ! ^°' r , where fix is an arbitrary unit vec- 
tor. Assuming a slow variation of staggered magnetiza- 
tions in the low-period spiral, we parameterize the or- 
der parameters as: L;(r) = e lqr </>i(r) fij(r), where rij(r) 
and 4>i (r) are the directions and magnitudes of the stag- 
gered magnetizations. The hard constraints |Sj| = S 
require that three unit vectors are orthogonal to each 
other and <j>\ + </>| + <j>\ = l. 27 Note that the phase factor 
e «Qo r _ _|_^ changes sign between consecutive layers of 
tetrahedra along the z direction (cf. the two layers of 
green tetrahedra in Fig. 3). 

The calculation can be further simplified in the J — > 00 
limit. The vanishing of the total magnetization on tetra- 
hedra in this limit requires 63 || d y hi, and the 4> fields 
can be analytically expressed in terms of the director 
fields hi. The effective energy functional in the gradient 
approximation can be solely expressed in terms of the 
director fields fix and h.2 \\ d y n.i x fix 27 . It contains two 
contributions: T = J dr (Ek + Edm)- The first term, 
originating from magnetoelastic coupling, gives penalties 
to variations of the staggered magnetizations: 



E K = 



K u 



(a c nx) 2 + (<9. y n 1 ) 2 + 2 («9 2 fix) 2 - (n 2 • 3A) 2 



(3) 

while the DM term contains Lifshitz invariants which are 
first-order in the gradients of fix : 



Dhi ■ (ax c^fix +b x d y hi — cx <9 z fix) /4. (4) 



Here we orient the axes of coordinate reference frame x, y, 
and z along the principal axes a, b, and c of the crystal. 

The form of the DM functional suggests spiral solu- 
tions with fix rotating about one of the principal axes and 
staying in the plane perpendicular to it. These special 
solutions were first obtained in Ref. 27. For example, the 



solution fix = (cos#(y), 0, smd(y)) with 8(y) = 27rSy + ip 
describes a helical order with coplanar spins lying in the 
ac plane (see Fig. 3). This spiral magnetic order pro- 
ducing a Bragg scattering at Q = 2w(0, 5, 1) is consistent 
with the experimental characterization of CdCr204. 3 

The energy functional T actually possesses a continu- 
ous symmetry related to the 0(2) rotational invariance of 
the helical axis. To describe the general coplanar spiral 
solution, we first introduce two orthogonal unit vectors 
lying in the xy plane: 

f = + cos $ x + sin $ y, 

C = — sin $ x + cos $ y (5) 

A coplanar spiral with spins rotating about the £ direc- 
tion has the form 

ni(r) = cos <?(£) C + sin 0(f) z, f = r ■ £ (6) 

Substituting ni(r) into the Eqs. (3) and (4) yields an en- 
ergy density: E[0] = —Dd^O + \K U (d^9) 2 , whose mini- 
mization gives a helical wavenumber 

Q = 2tt5 = d8/d£ = 2D/K U . (7) 

The 0(2) invariance indicates that the magnetic energy 
is independent of the angle 3>, i.e. (^[ni]/^ = 0. It is 
worth noting that this degeneracy is similar to the 0(3) 
symmetry for the helical axis observed in the spiral order 
of MnSi. 31 More importantly, as will be discussed below, 
it is exactly this additional symmetry that gives rise to 
the anisotropic dispersion of helimagnons. 

The coplanar spirals described above (with its axis 
lying in the xy plane) are also accidentally degenerate 
with a spiral in which spins rotate about the z axis. 27 
In real compounds, this accidental degeneracy as well 
as the aforementioned 0(2) symmetry are lifted by fur- 
ther neighbor interactions and cubic symmetry field of 
the lattice. In particular, the further-neighbor exchange 
contributes to a gradient term ( J 3 — J 2 )(<9 z rii) 2 in the gra- 
dient expansion. 27 The experimentally observed (0,6,1) 
spiral is selected by a large antiferromagnetic third- 
neighbor exchange J3 which favors coplanar spirals rotat- 
ing about the x or y axes. The explicit magnetizations 
of a spiral with the pitch vector parallel to the y-axis are 

S M (r) = a sin(Q ■ r + ip^) + c cos(Q ■ r + tp M ), (8) 

where \l = 0, 1, 2, 3 is the sublatticc index, tpo = ip, ifi = 
ttS + ip, if2 = it + ip, and ip^ = tt(1 + 5) + ip, with ip an 
arbitrary constant related to the U(l) symmetry of the 
coplanar spiral, i.e. 6(r) = Q ■ r + ip. 

It is worth noting that the pitch of magnetic spirals 
induced by DM interaction is usually rather long with 
the wavenumber of the order S ~ D/J ~ 1CP 3 — 10~ 2 . 
The relatively short pitch of spirals in CdCi'204 is due 
to the ineffectualness of J in stabilizing the magneitc 
order. Instead, the DM interaction competes with the 
magnctoclastic coupling and the helical wavenumber is 
given by S — D/K ps 0.1. 
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FIG. 3: (Color online) Coplanar helical order with spins lying 
in the ac plane and rotating about the b axis. The magnetic 
order corresponds to the solution fii = (cos 6(y), 0, sin 6(y)), 
with 0(y) = 2nSy + const. It produces a Bragg peak Q = 
2ir(0,5, 1) in neutron scattering, consistent with the experi- 
mental measurement on CdCr2 04. 

B. Spherical approximation 

While the above continuum theory gives a rather con- 
sistent and simple picture of helical magnetic orders in 
tetragonal pyrochlore antiferromagnets, it is based on the 
assumption of an infinite J. To go beyond this limit and 
investigate the magnetic orders when the various param- 
eters are of similar order, we resort to the spherical ap- 
proximation, in which the local length constraints |S j| = 

S are replaced by a global one J2i=i |Si| 2 = NS 2 , where 
N is the number of lattice sites. With this soft constraint, 
the model Hamiltonian (1) can be diagonalized with the 
aid of Fourier transform Si = ^ k S M (k) e lk ' r . Here we 
label a lattice site as i = (fi, r), where \i = 0, 1, 2, 3 is the 
sublattice index and r denotes its position. The Hamil- 
tonian then becomes 

n = j;x;^(k)s M (k).s,(-k) 
+ EEE a#(k)s-(k)s£(-k), (9) 

where the matrices J MI ,(k) and D^(k) are the Fourier 
transform of the exchange and DM interactions, respec- 
tively. The explicit form of the matrices are given in 
the Appendix A. The ground-state energy is given by 
the minimum eigenvalue A m ; n of the matrix — 
Jfj,L>3 a ^ + . In particular, the ordering wavevector 
Q is obtained by minimizing A m ; n (k) with respect to k, 
i.e. E = iVA min (Q). 

Fig. 4 shows the ordering wavevector Q as a function 
of DM interaction and third- neighbor coupling J3. We 
find three different helical orders characterized by an in- 
commensurate wavenumber 5 ~ D/ K u . At small J3, the 
ordering wavevector Q = 2tt(Q, 0, 1 — 8) indicates a mag- 
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FIG. 4: (Color online) Phase diagram of helical magnetic or- 
ders on pyrochlore lattice based on spherical approximation. 
Both D and J3 are measured with respect to the NN ex- 
change J, and we have set the odd-parity exchange anisotropy 
K u = 0.1J and 0.2 J while leaving the even-parity exchange 
anisotropy K g = 0. The green, yellow, and white regions 
are characterized by a ordering wavevector Q = 27r(0, 5, 1), 
27r(0, 0, 1 — 5), and 2tv(S, S, 1), respectively. 




FIG. 5: (Color online) A schematic diagram showing the out- 
of-plane corrections to a coplanar helical order. The coplanar 
spins are described by a vector field fii = (cos 6(y), 0, sin 6(y)), 
with 9(y) = 2n5y + const, while the out-of-plane components 
are given by Eq. (13). The inset shows the out-of-plane spin 
component AS as a function of r computed for the following 
set of parameters: J = 1.35 meV, K u = 0.21 meV, K g = 
meV, D = 0.14 meV, and J 3 = 0.28 meV. 



netic spiral in which spins rotate about the z axis. This 
spiral solution has been discussed in Ref. 27 and is acci- 
dentally degenerate with the coplanar spirals described in 
Eq. (8) in the J — > 00 limit. A new spiral order with or- 
dering wavevector Q = 2tt(6, 5,1) is obtained for a small 
window of J 3 (depending on the value of K u ). This spi- 
ral state with its axis pointing along the diagonal in the 
xy plane is also included in the general coplanar spiral 
solution Eq. (6) in the infinite J limit. As J3 is further in- 
creased, the (0, 5, 1) spiral observed in CdCr 2 04 becomes 
the ground state. For a helical wavenumber <5 ~ 0.1, this 
spiral order is indeed the magnetic ground state for a 
J3 rj 0.3 J, consistent with the ab initio calculations. 27 ' 34 
Interestingly, for small exchange anisotropics, e.g. 
K u = 0.1J, the spiral order (0,5, 1) with axis along ei- 
ther a or b axes is sandwiched by the spiral state with 
axis along the [110] direction, and the ordering wavevec- 
tor Q switches back to 2w(5, 5, 1) for largest J3 [Fig. 4(a)]. 
On the other hand, for larger anisotropy K u , the ground 
state remains the (0,5, 1) spiral for large J3. 

C. Perturbation correction to coplanar spirals 

We now focus on the helical magnetic order with 
wavevector Q = 2tt(0,5,1), which is relevant for the 
case of CdCr204, and use perturbation method to obtain 
corrections s, ~ (D(D/J) to the coplanar spiral solution 
Eq. (8) when the NN exchange is finite. In the ground 
state, the magnetic moment is parallel or antiparallel to 
the effective exchange field 

H» = dn/dS, = (Jij Sj - D« x Sj) , (10) 

3 

such that the torque T< = Sj x H; vanishes. Math- 
ematically, our perturbative calculation is based on a 



1/J expansion of the equilibrium equation T» = 0. The 
zeroth-order result corresponds to the J — > 00 limit and 
is given by Eq. (8). We express the exchange field as 
H; = Hi + the condition of zero torque requires 

(S< + Si) x (H, + = 0. 

Collecting only terms linear in D/J from this equation, 
we obtain a set of coupled equations: 

^ ./,,(•„>(/,, ;< <! II. (11) 

3 

]T \n;, cos 0, - Dfj sin 0, - ./,.. [cos % S f -*?]} = o, 

3 

(12) 

where 0, = Q • r + ?/v , an d % = @i — Sj. Sine the NN ex- 
change J can be expressed as (J/2)|M(r)| 2 , the J — > 00 
limit imposes hard constraints M(r) = on all tetrahe- 
dra. A finite value of J is thus expected to introduce a 
finite total magnetization. For simplicity, here we con- 
sider anisotropics dominated by a tetragonal distortion 
with odd parity only. Solving the above equations yields 
( M = 0,---,3) 

fey T~J 

s„(r) = - -cos(Q-r)b. (13) 

M ' (4J + 10K U ) y ^ ' K 1 

Interestingly, the first-order corrections does not modify 
the spiral pitch; the incommensurate ordering wavevector 
is still given by the zeroth order expression Q = 2ir5. 
As the correction is uniform for the four sublattices, the 
finite value of J gives rise to a modulated ferromagnetic 
component M(r) on individual tetrahedron. 
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FIG. 6: (Color online) The energy decrease AE(t) = E(t) - 
E(t = 0) as a function of simulation time starting with (a) 
coplanar spiral solution Eq. (8) and (b) initial state including 
first order correction (13). The energy and time are mea- 
sured with respect to J and 1/J, respectively. We performed 
the simulations on a lattice with 16 x 10 3 spins with periodic 
boundary conditions. Panels (c) and (d) show the correspond- 
ing averaged spin deviations es(i) = J2i |St(<) — S>|/JV from 
the coplanar state and the inclusion of the first order correc- 
tion respectively as a function of simulation time. 



D. Relaxation dynamics simulations 

While the perturbation correction gives a better ap- 
proximation of the helical order in pyrochlore lattice, a 
high-precision description of the magnetic structure is 
required for the calculation of spin-wave excitations. A 
common approach to obtain the magnetic ground state 
is the low-temperature Monte Carlo simulations. How- 
ever, the spin configurations obtained from Monte Carlo 
simulations are still too noisy for the purpose of the cal- 
culation of magnon spectrum. Instead, here we resort 
to the overdamped dynamics simulations based on the 
Landau-Lifshitz-Gilbcrt (LLG) equations. In particular, 
since the coplanar spiral solution with the perturbation 
corrections already provide a close approximation to the 
true ground state, such relaxation simulations can effi- 
ciently bring the system to the desired magnetic order. 

The LLG equation is a first-order differential equation 
describing the dynamics of classical magnetic moments: 



= S. t x Hi 



S ' 



(14) 



O 1-0 




where H^ is the effective exchange field defined in 



FIG. 7: (Color online) Structure factor of the helical spin or- 
der obtained from the relaxation dynamics simulation. Sharp 
peaks appear at the ordering wavevector Q = 27r(0, 5, 1). 



Eq. (10), and Qg is a dimcnsionlcss damping parame- 
ter. We employ a finite-difference method to integrate 
the LLG equation numerically. 35 The technical details of 
the numerical method can be found in Appendix B. We 
performed our simulations on a cluster of 10 3 cubic unit 
cells with TV = 16000 spins. Periodic boundary condi- 
tions were used in all simulations. We have also tested 
our algorithms on some well studied Hamiltonians and 
obtained excellent agreement with exact solutions. 

We start our dynamics simulation from the coplanar 
spiral state both with and without the first-order correc- 
tions. Fig. 6 shows the simulation results for parameters 
K u = 0.1J, D = -0.069J and J 3 = 0.2J; the corre- 
sponding helical wavenumbcr is 5 = 0.1. In both cases, 
the energy monotonically decreases with simulation time, 
indicating a steady relaxation towards the true helical 
magnetic ground state. Also shown are the averaged de- 
viations of spin configurations Sj(t) from the coplanar 
helical state (with and without first-order correction) as 
a function of simulation time. The rather small devia- 
tions of the order 10~ 2 indicate that both the coplanar 
spiral and the first-order result are indeed already good 
approximations to the ground state. Fig. 7 shows the 
structure factor 5(k) = |(1/A) J^i Si exp(ik • r^) | 2 of the 
spin configurations obtained from the relaxation dynam- 
ics simulations. As expected, the structure factor ex- 
hibits peaks at the ordering wavevector Q = 27r(0, S, 1). 



III. SPIN- WAVE EXCITATIONS 

In this Section we study magnetic excitations of the 
helical magnetic order on the pyrochlore lattice. We first 
discuss the so-called helimagnons in the continuum limit. 
The appearance of this special gapless mode originates 
from the unique symmetry-breaking in chiral magnets. 
However, since the spiral pitch in CdCr204 is relatively 
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short compared with helical orders in other systems, the 
crystal effects are expected to modify the helimagnon dis- 
persion significantly. We thus also compute the spin-wave 
dispersions based on microscopic spin models. Starting 
from the high-precision ground state obtained in the pre- 
vious section, we performed exact diagonalizations on 
a large unit cell to study the nature of the low-energy 
acoustic-like modes. The results are then compared with 
the predictions based on helimagnon theory. We also em- 
ployed dynamics simulations on hnite systems to obtain 
the full magnon spectrum and compare the results with 
the experiments. 



A. Helimagnons in the continuum limit 

We first consider magnetic fluctuations in the helical 
order within the framework of the continuum theory. For 
a spiral with a fixed axis, the obvious soft modes are 
phase fluctuations associated with the coplanar director 
field. More specifically, let us consider a spiral with its 
axis along the y-direction: rii = (cos 8, 0, sin 8). Intro- 
ducing a slowly- varying perturbation to the phase angle 
8 = Qy + 4>{y)i the energy functional associated with 
the fluctuations is £[ip] = K u J dr [V"0(r)] 2 . This indeed 
describes a soft mode similar to that of a planar mag- 
net. However, as is well known in chiral nematics and 
smectic liquid crystals, a phase fluctuation of the form 
ip = a £ corresponds to a simple rotation of the helical 
axis. 29 Since coplanar spirals with different axes are en- 
ergetically degenerate (see discussions in Sec. II A), such 
rotations cannot cost energy, yet £ ~ a 2 ^ for this 
particular phase fluctuation. Consequently, there cannot 
be any (d^) 2 term in the effective energy density. The 
correct energy functional thus has the form 29 ' 31 



total energy functional becomes 



ro 



T[m, ip\ = j drm 2 (r) + £[ip] 



(17) 



Here ro is an effective mass of the zero-wavevector mode. 
The magnetization dynamics is governed by a generalized 
Landau-Lifshitz equation: 



S t S /1 (r) = - 7 S M (r) x 



ST 



(18) 



where 7 is a constant. Substituting (16) into the above 
equation and following similar steps outlined in Ref. 31, 
we arrive at the following coupled equations 



d t m(r) = —7 



ST 



d t ip(r) = 7 



ST 
5m(r) 



(19) 



To obtain the spectrum of the soft mode, we eliminate 
m in the above equations to obtain a wave equation for 
the ip field 

d t 2 V(r) = -r o7 2 [~c z dl ~ c^ 2 + ^d^j V(r). (20) 

Assuming a space-time dependence tp ~ exp(ik-r — ie^t), 
the energy of the soft mode is 



(21) 



The dispersion is highly anisotropic: for wavevectors par- 
allel to the spiral £-axis direction or the z-axis the dis- 
persion is linear as in an antifcrromagnet, while it is 
quadratic for wavevectors parallel to the x-axis. 



£ = -ldr 



(15) 



where cn, c± ~ K u , c z ~ J3 are elastic constants, and 
Q = 2ttS is the pitch wavenumber. Note that the first 
term in the above equation arises from a rather large J3 
which gives a penalty to variations along the z direction. 

A derivation of the dynamics equations starting from 
the microscopic model (1) is tedious and almost in- 
tractable. To obtain the dispersion of the soft mode, it is 
more illuminating to employ a simple phenomenological 
approach based on the time-dependent Ginzburg-Landau 
theory 29,31 In this formalism, the Goldstone mode ip cou- 
ples to a mode m at zero wavevector, which is soft due 
to spin conservation. Schematically the coarse-grained 
magnetization field is 

S„(r) ~ S M (r) + m(r)y (16) 
+ tp(r) [cos(Q ■ r + ip^) x - sin(Q • r + tp^) z] , 

where S M is the coplanar spiral given in Eq. (8) and ip^ 
denotes the relative phase shift of the sublatticcs. The 



B. Exact diagonalization with large unit cell 

As discussed in Sec. II A, the helical wavenumber in 
tetragonally distorted pyrochlore is determined by the ra- 
tio of DM interaction to exchange anisotropy (instead of 
the NN exchange constant). The spiral pitch in CdC^CU 
is relatively short compared with that in other helical or- 
ders. We thus expect significant anisotropy coming from 
the microscopic crystal fields. In addition, the (0,(5,1) 
spiral which is relevant to our case is stable only in the 
presence of a small third-neighbor coupling. A gradient- 
expansion analysis shows that spin stiffness along the z 
axis is increased in the presence of J3. To address these 
issues, here we perform exact diagonalizations with a 
large unit cell to investigate the low-energy acoustic-like 
magnons and compare the results with predictions based 
on the helimagnon theory. 

To begin with, we first choose a commensurate unit 
cell composed of A cubic unit cells along the b axis (each 
cube contains 4 tetrahedra), corresponding to an order- 
ing wavevector Q = 2n(0,8, 1) with 6 = 1/A. The DM 
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FIG. 8: (Color online) Spin-wave dispersions along three principal axes [panels (a)-(c)] obtained from numerical exact diago- 
nalization of a large unit cell. The results are fitted to the phenomenological dispersion e k = \J ctk 2 + /3k 4 . The magnon energy 
is measured in unit of JS. Panels (d)-(f) show the dispersion parameters a and /3 as a function of the helical wavenumber S. 
These data are fitted to a =const, a oc 8 2 , and f3 cx S~ 2 , respectively. 



constant that is required to stabilize such a spiral is then 
determined using the expression 

D = -(3/V2)K u t&mrS (22) 

This is the discrete version of Eq. (7) for a given fixed 
anisotropy K u , and is obtained by substituting the copla- 
nar spiral ansatz (8) into Hamiltonian (1) and minimizing 
the resulting expression with respect to 5. The relaxation 
dynamics simulation discussed in Sec. II D is carried out 
to obtain an accurate description of the ground state. We 
then introduce small deviations to spins in the ground 
state: S, = S? + <x,, where (Ji _L S° denotes transverse 
deviations. Note that because of this constraint, there 
are two degrees of freedom associated with <x^ at each 
site. The LLG equation is linearized with respect to the 
small deviations erf. 

^ S, 1 ;>>/, ; <t. I),, <t ; i II 1 ; a,. (23) 

3 

where Jij includes J, Kij and further-neighbor ex- 
changes, H° = Ylj(JijS < j — x S°) is the equilibrium 
local exchange field, and we have set the damping con- 
stant olq to zero. The eigenmodes of the LLG equation 
has the form cr, = cr m cxp(ik ■ r — ie-^t), where the site 
index i = (r,m), where r denotes the position of the ex- 
tended unit cell containing A cubes and m labels spins 



within this block. Note that the above eigenmode {<r m } 
has a momentum Q + k due to the underlying struc- 
ture of the supercell. Arranging cr m = {a'^ n ,af n ) into 
a column vector u = (erf , a\, ■ ■ ■ , a* , erf) of dimension 
?i = 2 x 4A x 4, the linearized LLG equation is recast 
into an eigenvalue problem T(k) • u = — ie k u. Numeri- 
cal diagonailization of the n x n matrix T(k) gives the 
dispersion of the low-energy magnons in the vicinity of 
ordering wavevector Q. 

The resulting magnon dispersions along the three prin- 
cipal directions are shown in Fig. 8(a)-(c). The number 
of spins in the unit cell is defined by the value of the 
spiral pitch 5. We have considered 6 = 0.2, 0.1, 0.05, 
which gives N s = 80, 160, 320, respectively. The numer- 
ical spectra are fitted to a phenomenological dispersion: 
e k = sj ah 2 + /3fc 4 . The dispersion parameters a and (3 
obtained from the fitting are shown in Fig. 8(d)-(f) as 
a function of the helical wavenumber S. First, we note 
that the dispersion along the y (helical direction) and z 
(staggering direction) axes can be well approximated by a 
linear relation, e k ~ vk, with v = \fa weakly dependent 
on S, see Fig. 8(d). The large spin- wave velocity along 
the z axis can be attributed to a large J3 , consistent with 
the analysis of gradient expansion. 2 

On the other hand, the spectrum along the x-axis ex- 
hibits a predominant quadratic behavior at small fc, sim- 
ilar to that of a ferromagnct. Although the continuum 
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theory predicts an exact quadratic dispersion along the 
x direction, Eq.(21), careful fitting shows a small linear 
component a, which can be attributed to the discrete 
cubic symmetry of the underlying lattice model. Di- 
mensional analysis indicates that a should scale as Q 2 , 31 
while the coefficient (3 ~ 1/Q 2 according to the heli- 
magnon dispersion Eq. (21). These scaling relations are 
indeed confirmed by our exact diagonalization as demon- 
strated in Figs. 8(e) and (f). 



IV. MAGNON SPECTRUM OF CdCr 2 4 
A. Linear dynamics simulation 

Although the exact diagonalization method discussed 
above provides a rather accurate magnon spectrum, the 
complicated band-folding due to large unit cell makes it 
difficult to compare the results with experiment. Instead, 
here we employ a less accurate but more direct approach 
based on simulating the real-time dynamics of spins in a 
large finite lattice. Specifically we introduce a local basis 
for the transverse fluctuations: er^ = af ef +a\ e*, where 
ef are two unit vectors orthogonal to the equilibrium spin 
direction . Substituting <Tj(t) into the linearized LLG 
equation yields 

dtof = (rf - a G t?) /{I + a%) + Cf (t), 

(24) 

d t a\ = {a G rf + T ?)/(l + a%) + (?(t), 

where rf and are linearized torques projected onto the 
local basis; see Appendix C for details of the calculations. 
A small damping is added to ensure stable simulations. 
The dynamics simulation is initiated by a short pulse 
d(t) = ( n cxp(—t 2 /w 2 )Sifi localized at r 2 = 0. Here Co 
determines the strength of the initial perturbation and w 
is the width of the pulse. 

The differential equations (24) are integrated with the 
aid of the fourth-order Runge-Kutta method. The solu- 
tion (Ti = <r m (r, t) describes the evolution of each indi- 
vidual spin due to the localized short pulse. Since the 
system is driven by a white source with flat spectrum, 
we expect magnetic excitations of various energy and 
momentum are generated in our simulations. The spin- 
wave spectrum is then obtained via the Fourier trans- 
form of the simulation data. First we performed the 
spatial Fourier transform <x m (k, t) = 5Z r o" m (r, t)e zk ' T , 
where only those momenta k which satisfy the periodic 
boundary conditions of the finite lattice are allowed. In 
Fig. 9 we present the time evolution of cr m (k, t) for dif- 
ferent values of k. The magnon spectrum is then ob- 
tained the temporal Fourier transformation: <r m (k, e) = 
/ e le *<r m (k, t)dt. The spectrum obtained for a set of 
model parameters that are realistic for CdCr204 is dis- 
cussed in the next subsection. 
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FIG. 9: (Color online) Time evolution of uq (k, t) for different 
values of k. 

B. Comparison with experiments 

The magnon spectrum of CdCr204 has been studied 
by Chung et al. using inelastic neutron scattering. 3 Here 
we present our numerical calculations and compare the 
results with the experimental findings. The fact that 
the helical wavenumber S « 0.09 can be used to fix the 
ratio of the DM interaction to the exchange anisotropy 
K u . On the other hand, values of K u and K g can be 
estimated by the magnon energies at the zone center, and 
we find that even tetragonal distortion does not play a 
significant role here. As for further-neighbor exchanges, 
ab initio calculations find a negligible J2 and a quite large 
J3 with antiferromagnetic sign. 27,34 Using the following 
set of parameters: J = 1.46 meV, K u = 0.34 meV, K g = 
meV, D = 0.078 meV, and J3 = 0.30 meV, we obtained 
very good agreement with the experimentally measured 
spectrum, as shown in Fig. 10. 

V. CONCLUSION AND DISCUSSION 

We have studied the spin-wave excitations in a spe- 
cial helical magnetic order on the pyrochlore lattice. 
This spiral structure characterized by a wavevector Q = 
27r(0, 6, 1) is observed in frustrated spinel CdCr204. Mo- 
tivated by the experiments, we assumed that magnetic 
frustration is completely relieved by a tetragonal dis- 
tortion which preserves the lattice translational symme- 
try. We have employed various analytical approaches 
and numerical methods to investigate the interplay of 
Dzyaloshinskii-Moriya interaction and further-neighbor 
couplings. Our spin- wave calculation with parameters 
inferred from structural data and ab initio calculations 
agrees well with the measured magnon spectrum. Com- 
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FIG. 10: (Color online) (a) Spin wave spectrum of a Neel 
state with wavevector q = 2n(0, 0, 1). This collinear order is 
the ground state of H with parameters J = 1.35 meV, K u = 
0.21 meV, D = meV, and J 3 = 0.28 meV. (b) Fitted spin 
wave spectrum to the neutron scattering measurement (white 
circles and rectangles) 3 . The final set of parameters are J = 
1.35 meV, K u = 0.21 meV, D = 0.14 meV, and J 3 = 0.28 
meV. Note that adding DM interaction lifts the gap at q = 
27r(l, 1, 0) and splits each band into three. All the parameters 
are determined based on this fitting. Note that vertical peaks 
at the ordering wave vector q = 27r(l, 1 ± S, 0) are numerical 
artifacts. 



bined with a systematic analysis of the helical structure, 
our study helps clarify the underlying mechanisms that 
stabilize the observed helical order and provides a quan- 



titative description of the model Hamiltonian. 

We have also shown that the continuum model of the 
spin spirals based on a gradient-expansion of the or- 
der parameter exhibits a continuous symmetry similar to 
that of smectic liquid crystals. A Goldstone mode with 
highly anisotropic dispersion results from the unique bro- 
ken symmetry of the helical order. Although the emer- 
gent continuous symmetry is broken by the crystal fields 
effects at the microscopic level, our exact diagonalization 
on a large unit cell finds acoustic-like magnons which are 
reminiscent of the Goldstone mode discussed above. 

Several recent experiments have established CdCr204 
as the paradigmatic example of magnetic ordering in- 
duced by spin-lattice coupling in geometrically frustrated 
magnets. Although an overall tetragonal distortion with 
a = b < c is observed below the ordering temperature, 
details of the lattice distortion remain to be clarified ex- 
perimentally. The issue at hand is whether the distor- 
tion breaks the lattice inversion symmetry. As pointed 
out in a previous theoretical study, 27 the broken parity 
plays an important role in stabilizing the observed he- 
lical order. This is because the lack of inversion sym- 
metry endows the crystal structure with chirality, which 
is transferred to the magnetic order through relativistic 
spin-orbit interaction. Experimental evidence implying 
a broken parity in CdCr2 04 was recently reported in in- 
frared absorption spectra. 14 Our numerical calculation 
on spin-wave dispersions also indicates a predominantly 
tetragonal distortion which breaks the inversion symme- 
try, giving further support to the scenario that the helical 
magnetic order inherits its chirality from a lattice with 
broken parity. 
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Appendix A: Spherical approximation 

The matrices J M „(k) and D"^(k) in the spherical approximation (Eq. 9) are defined as follows: 



J(k) 



/ ^ J 3 (C X y + C XZ ~ + Cy Z ) (2 J - K u 



3A u Sy £ (2,7 Ky^jCxz ~f~ 3K u s X z 



2(</ ~\~ K u ^C X y 



\ 



(2 J - K u )cy Z - 3K u s yz 2J 3 (c xy + c xz + c yz ) 2(J + K u )c xy (2 J - K u )c xz + 3K u s xz 

(2t7 I^u)c X z 3K u S X z 2(t/ -)- Ku)Cxy 2^/3 (c X y -\- C xz ~ Cyz*) (2*7 Ku)Cyz 3K u Syz 

\ 2(t7 + K u )c X y (2J K u )c xz 3K u s xz (2J K u )c yz - 3K u Syg 2J3 (c X y + Cyg + c xz ) ] 



(Al) 
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kp), c 


af) = 


cos(fc a 



(A2) 



cos(2(fc Q + fc^)), and s a( 3 = isin(k a + hp). 

Appendix B: Finite-difference method for integrating LLG equation 



In this appendix, we discuss a finite difference scheme for numerical integration of the Landau-Lifshitz-Gilbert(LLG) 
equation by adopting the method introduced by Serpico et aP 5 . The LLG equation can be written in the following 
form: 



dS, ( ct G dS t 

Note that the effective exchange field H; is defined in the Eq. (10). 

Then the following approximations are applied on the midpoint between time step n and n + 1 



(Bl) 



-Qj-(t n +l/2) - 



Si (t n ) 



At 



0(At 2 ) 



(B2) 



(B3) 



where At = t n+ i — t„ and t n+ i/ 2 denotes midpoint between t„ and 

We now express EL(t„ +1/ / 2 ) in terms of Hj(i n ) and Hj(f n _i) by applying the midpoint approximation (B3) twice : 



3 1 

Hi(t„ + i/ 2 ) = -Hj(t„) - -H,(t n _i) 

By substituting Eq. (B2),(B3), and (B4) into Eq. (Bl) we get: 



(B4) 



Si(t n +l) — S;(t„) 
At 



Sj(t n +i) + Si(t r; 



-H,(t„) - -Hi(t„_i) + — — 



(B5) 



which yields the expression for Sj(t„_|_x) m terms of previous time step spin components S,-(t„) and Sj(t n _i): 



Si(t n +i) — Sj(t n +i) x 



hli{t n ) - ^Hiitn-x)) At - |s<(t„) 



(B6) 



( jHi(t„) - ^Hi^x) 



(B7) 
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Then the spin configuration Si (t n +i) is obtained by solving the matrix equation: 



Si(tra+l) — Bj(t n , t n — \) Aj (t n , t n — \) Si (t n ) i 

where Aj(t„, f„_i) and (t n ,t n —i) are both 3x3 matrix defined as follows : 



(B8) 



i -/c? K.y 

fC? 1 -£? 
-AC? Kf 1 



(B9) 



Bi(t n ,t n _i) 



i n -n 

-T? 1 J? 
T? -Tf 1 



where & = [%Hi(t n -i) - |H 4 (t„)] At and P t = t t + ^Si(t„). 



(BIO) 



Appendix C: Torques in the linearized LLG equation 

The torques in the linearized time-dependent LLG equation (24) are calculated in this Appendix. First we choose 
the following paramctrization for the local basis: 



ef = cos 9i cos 4>iii+ cos 6^ sin cj>ih — sin 9iC , 

ej = — sin <j^a + cos </>;b , 

fii = sin 0j cos + sin 0j sin <pih + cos 0jC , 



(CI) 



where Qi and ^ are polar and azimuthal angles that determine the orientation of spin at site i. Introducing the 
following notations for simplicity: Cj = cos(9j, Sj = sin#.j, Cj = cos0j, s, = sin<^. the linearized torques projected 
onto the localized basis are: 



3 3 



(C2) 



3 

X]|( /) ' / ; r ' ~ D ij S i) C * S 3 ~ ( D ij S i + D fj^) r J r J + ( D ij c iCi + D ij s i) r , s J + J -J {<kCjCi-j + SiSj) 



^2[ D ij C iCi-j + (DfjCj + DVjSj) S. 



(C3) 
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